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ABSTRACT 

Wave travel-time shifts in the vicinity of sunspots are typically interpreted as arising pre- 
dominantly from magnetic fields, flows, and local changes in sound speed. We show here that 
the suppression of granulation related wave sources in a sunspot can also contribute signifi- 
cantly to these travel-time shifts, and in some cases, an asymmetry between in and outgoing 
wave travel times. The tight connection between the physical interpretation of travel times and 
source-distribution homogeneity is confirmed. Statistically significant travel-time shifts are re- 
covered upon numerically simulating wave propagation in the presence of a localized decrease in 
source strength. We also demonstrate that these time shifts are relatively sensitive to the modal 
damping rates; thus we are only able to place bounds on the magnitude of this effect. We see 
a systematic reduction of 10-15 seconds in p-mode mean travel times at short distances (~ 6.2 
Mm) that could be misinterpreted as arising from a shallow (thickness of 1.5 Mm) increase (~ 
4%) in the sound speed. At larger travel distances (~ 24 Mm) a 6-13 s difference between the 
ingoing and outgoing wave travel times is observed; this could mistakenly be interpreted as being 
caused by flows. 

Subject headings: Sun: helioseismology — Sun: interior — Sun: oscillations — waves — hydrodynamics 

1. INTRODUCTION 



The discovery that sunspots support oscillations (e.g., see the review by iBogdan fe Judgjl2006h was 
important for it introduced the possibility of using measurements of phase-shifts in the propagating waves to 
probe the underlying structure and dynamics of these enigmatic objects. Some of our current observational 
understanding o f the sunspot inter i or comes from invers e theory applied in conjunction with time-distance 



helioseismology ( Duvall et al. 19931 : Gizon fc Birch 20051) on waves in these regions. S ubsequent to the stud- 
ies of flows in and around sunspots bylDuvall et a" . (1996). inversion s utilizing the ray |Kosovichev fc Duvalll 
1997 ). Rytov ( Jensen fc Pijpers 20031 ). a nd Born (Birch et al. 2004) approximations were performed to re - 



cover the interior structure of sunspots ( Kosovichev et al. 2000 : Jensen et al. 2003t Couvidat et al. 2006f h 



In recent years, many of these results have come into question because the analyses do not explicitly ac- 
count for the influence of magnetic fields on wave travel times. Apart from direct mechanical effects on the 
waves, magnetic fields are also responsible for impeding the action of near-surface convection, commonly 
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beli eved to be the source of waves (e.g., IStein fc Nordlundl 12000 ) . Despite the work of IWoodardl (|1997h 
and lGizon fc Birchl (|2002f) . causal factors of travel-time shifts such as wave damping and source distribution 
inhomogeneity in the context of sunspots have not been studied in detail. 

Gizon fc Birchl |2002h first derived the linear sensitivity of /-mode travel ti mes to local changes in source 



strength, later corroborated through time-distance analyses of artificial data bv lHanasoge et alj (|2007| ). The 
concept of variations in source-strength engendering travel-time shifts can be somewhat mystifying. Surely 
waves do not speed up or slow down when a source emits a wave of half the amplitude, as the naive 
interpretation seems to indicate? The answer lies in the manner in which travel times are computed; stripped 
of physical interpretation, travel times are obtained by fitting cross correlations of velocity (or intensity) 
signals between pairs of points or a point and an annulus. The measurement points do not constitute a 
source-receiver pair as in the typical geophysical situation; rather, all waves that contain coherent phase 
information at these points contribute to the cross correlations. The wave travel times measured in a system 
with a spatially uniform distribution of sources and a specific set of damping rates have certain expectation 
values. However, it is conceivable that over a region where the directionality or spatial distribution of waves is 
biased, the contributions by wave packets (to the cross correlations) from disparate directions and points are 
not in the same proportion as in the spatially uniform case. Consequently, there is a shift in the expectation 
value of the travel time in this region. In fact, the term 'travel time' is better interpreted as a quantity that 
describes the statistics of the wave field than the physical wave travel time between the measurement points. 
Damping also plays an important role, for it determines the extent of coherence of the waves and the degree 
of contribution to the cross correlations. These can be serious issues in sunspots, because of the possible 
lack of sources and the putative excesses in damping and absorption (e.g.. lBraun et al.lll987l ). 



Mean travel times are defined as the average of the in- and outgoing wave travel times, while difference 
travel times are obtained by subtracting the two. We posit that the classical interpretation of mean travel- 
time shifts as mostly arising from changes in the sound speed and difference travel-time shifts predominantly 
from flows in sunspots is incomplete because the lack of wave sources can also cause significant mean and 
difference travel-time shifts; this effect is demonstrated here via numerical simulations and semi-analytical 
methods. In ^ we describe the numerical machinery employed to perform the simulations and discuss 
the impact of horizontal boundary conditions on the resultant time s hifts. In order to ch aracterize the 
influence of damping rates, we apply the semi-analytical techniques of Gizon fc Birchl (2002). We analyze 
the simulated data with methods of time-distance helioseismology in comparisons are drawn between the 
results of simulations and the semi-analytical models. Finally, we summarize and conclude in §3J 



NUMERICAL PROCEDURE AND TEST CASES 



Using techniques developed in lllanasoge et al] ( 2006h . Hanasoge et al. ( 2007 ). and Hanasoge ( 2007a ). 
wave propagation in the near-surface layers of the Sun is simulated in a box of dimension 400 x 400 x 35 Mm 3 , 
where th e third dimension is depth. The background stratification is a convectively stabilized form of 



model S (jChristensen-Dalsgaard et al.l I1996T ) , described in Appendix [X] Waves are stochastically excited 



by introducing a forcing term in the vertical momentum equation; the forcing function is prescribed such 
that a solar-like power spectral distributio n is obtained. T he solution is temporally evolved using a second- 
order optimized Runge-Kutta integrator (|Hu et al.l ll996). The vertical deri vative is resolved using sixth - 
order compact finite differences with fifth-order accurate boundary conditions (jHurlburt fc Rucklid gc 2000). 
Depending on the choice of boundary conditions, the derivatives in the horizontal directions are computed 
cither using these compact finite differences (absorbing conditions) or the Fast Fourier Transform (periodic 
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boundaries) . The validation and verification of the code is discussed in Appendix [B] 

The power spectrum and snapshots of the oscillation velocities derived from a 12 hour long 'quiet' 
simulation are displayed in Figures [1] and [2] respectively. To simulate source suppression, the forcing term 
is multiplied by a spatial function that mutes source activity in a circular region of 10 Mm radius (i.e., 
the forcing function assumes a reduced value in this region). Based on estimates of em itted energy flux in 
sunsp ot umbrae, which range from 10 - 20% of the average value in the quiet Sun (e.g.. ISchiissler fc Vogler 
2006), we perform two simulations, one with source strength in the disc region set to zero and another with 



20% of the 'quiet' value. The two simulations possess very similar travel-time maps; therefore, we only show 
results from the simulation where the sources were completely suppressed. 



Power spectrum 
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Fig. 1. — Power spectrum obtained from a 12 hour 'quiet' simulation. Some ridges have been labelled. The 
symbols mark independent estimates (obtained using MATLAB) of the eigenfrequencies of the altered solar 
model. The agreement between computation and theory appears reasonable. 

The power spectral distribution of oscillation modes and the steady-state wave flux emerge from an 
interplay between source activity, wave damping, and mode mass. The non-scattering nature of source- 
strength perturbations complicates matters because the phase measurements are sensitive to the degree of 
inhomogcneity, which in turn is dependent on the intensity of the ambient wave flux. One can imagine that 
in the limit of an arbitrarily large wave flux, the time-shift effects of the suppressed source may altogether 
vanish (or reach some asymptotically small value). It is therefore important to investigate this matter in 
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Fig. 2. — Snapshots of the normalized vertical component of the oscillation velocity {yfpoc v z ) - vertical and 
horizontal (at z — 200 km) cuts from a 'quiet' simulation (absorbing boundary case) are displayed. The 
units are arbitrary and same for both panels. Energy conservation requires an increase in velocities to offset 
the sharply decreasing density in the near-surface layers - and hence the choice of this normalization (po is 
the density and c the sound speed). 

some detail. The wave flux in the computations is set by the choice of boundary conditions and damping 
rates. All other parameters being identical, absorbing horizontal sides evidently engender a weaker ambient 
flux than their periodic counterparts; thus we may study the impact of boundary conditions on the time 
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shifts through numerical experiments with these choices of horizontal boundaries. In an indirect manner, 
these boundary conditions mimic higher and lower damping rates. In the case of the former, the lack of 
incoming waves from regions external to the boundaries sets the damping length (or maximum propagation 
distance) to approximately half the size of the computational region, since the perturbation is always at the 
center. This results in a dearth of p modes that travel large distances or those that possess long lifetimes. 
Contrarily, in the case of periodic boundaries, waves exit from one boundary only to re-enter the domain 
from another; if the modal damping rates are unrealistically small, these perpetually propagating acoustic 
waves will rapidly fill up the domain, thereby significantly diluting the effects arising from the suppressed 
sources. Roughly we may conclude that the low wave damping limit is given by the periodic case and high 
damping limit by the absorbing case. 

Unfortunately, due to poor observational constraints on damping rates, it is unclea r as to which of these 
situations is pre ferable. The linewidths recov ered from Michelson Doppler Imager (MDl lScherrer et al.lll995h 
observations by iKorzennik et al.l (|2004 ) and iDuvall et al.l (|1998l ) dif fer by almost a factor of 2. Moreover, 
the complicated functional dependence of damping on frequency ( e.g. IKorzennik et al.ll2004j ) makes it all but 
impossible to implement it in the computation. Thus we can only hope to place bounds on the extent of 
the effect of suppressed sources since precise estimates are closely tied to the non-trivial feat of accurately 
matching the simulated wave power spectral distribution with the observational one. The outcomes of these 
tests are discussed in 331 



2.1. Theoretical model 



In order to gain an appreciation for the effects of damping on the conclusions of this paper, we create 
semi-analytical forward models in the manner of Gizon & Birch (2002). These forward models predict the 
time shift associated with a specific perturbation. The starting point is the temporal Fourier transform of 
equation (22) from Gizon & Birch (2002), which gives the expected value of the cross-covariance, C, in terms 
of Green's functions Q and the source covariance M, 



C(n, r 2 , w) - (2ir) 2 J J ds £T (n, s, w)0*(r 2 , s, w)Af«(s, w) . 



(1) 



The integration variable s runs over the horizontal position of all the wave sources, r\ and r-i are the positions 
of the two observation positions, and lo is the temporal frequency. Notice that in writing this equation we 
have assumed that the wave sources are uncorrelate d in space. I n order to compute equation fl} we use the 



normal-mode summation of Green's functions from iBirch et al 
damping, one based on the line-widths measured by 



rates (approximately those measured by IDuvall et al 



( 200% which include two models of wave 



^orz ennik et al.l (12004) and the other with twice these 



1998). We use the source covariance from Birch et al 



(|2004l ). though modified to include the reduction of source strength inside the disc of radius 10 Mm. It is 
important to note that the type of source used in this particular forward model is quadrupolar in nature, 
whereas we employ dipoles in the simulation. With these ingredients, the expected value of the point-to- 
point cross correlation (Eq. [1]) can then be computed numerically and averaged to obtain center-to- annulus 
cross correlations. In we shall further discuss the connection between the horizontal boundary conditions 
implemented in the simulations and the damping rates incorporated in this theory. 
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TRAVEL TIMES AND POWER CORRECTION 



The p-mode travel times are measured using the procedure described in lCouvidat et al.l (J2006) . In order 
to estimate the travel times accurately, broad phase-spee d filters were implemen ted to avoid contaminating 
the first bounce ridge with the filter artifact (see Table 1 of lflanasoge et al~ ( 2007 ); the FWHM was multiplied 
by 4). The p-mode cross-correlation branches corresponding to positive and negative times (outgoing/ingoing 
waves) averaged over the source-perturbation area in comparison to the average cross correlation for the quiet 
simulation with absorbing horizontal sides are shown in Figure [3] for A = 24.35 Mm, where A = \ri — r%\ 
is the distance between measurement points. There is a noticeable asymmetry between the outgoing and 
ingoing waves, especially at larger distances where the outgoing waves contain almost all of the travel-time 
shift. Choosing the center of the source suppression to be the zero point, ingoing and outgoing travel-time 
shifts are azimuthally averaged and plotted in Figure 2J The reduction in the p-mode mean travel times 
seen in panel a of Figure 2] is comparable, magnitude- wise, to the 15 s increase (azim uthal average) seen fo r 
A = 6.2 Mm in some sunspots (NOAA 8243, from high-resolution MDI observations. I Couvidat et al.l l2006). 
The asymmetry between in and outgoing waves for the travel distance of A = 24.35 Mm results in significant 
shifts in the difference times, of the order of 12 s (panel b, Figured]). 



In contrast, the simulations with periodic boundary conditions show reduced shifts (also see lParchevskv et al 

2007), of the order of 10 s in the mean times for A = 6.2 Mm and 6 s in the difference times for A = 24.35 
Mm. Evidently, the reappearance of waves from the opposite boundary upon their exit from one has led 
to the prevalence of a larger wave flux in the computational domain. As pointed out in $2j the wave flux 
plays a crucial role in bounding the effect of non-scattering source perturbations. The situation is rendered 
interesting by the close correspondence between the theory of £12. II and the simulations, as seen between 
the upper and lower pairs of rows in Figure [4] Higher damping rates lead to larger time shifts and vice 
versa, analogous to si mulations with the abs orbi ng and period i c bou ndaries respectively. The conflicting 
linewidth estimates of iKorzennik et al.l (|2004h and lDuvall et alJ (j 19981 ) make it difficult to conclusively pick 
one simulation over the other. In fact, it is probably fair to say that realistic magnitudes of the time shifts 
lie somewhere between the estimates derived from the absorbing and periodic cases. 

If the simulations are believed to be representative of reality, and the travel times of in- and outgoing 
waves are appropriately 'corrected' to account for the possibility of missing wave sources in sunspots, we 
might expect a significant change in the mean travel times for A = 6.2 Mm. Moreover, the asymmetry 
between the in/outgoing waves seen for A = 24.35 Mm (ingoing ~ -10 s, outgoing ~ -40 s, azimuthal 
averages for sunspot NOAA 8243) could be reduced somewhat by applying these corrections. We show in 
Figure[5]that travel-time shifts associated with source suppression and sound-speed perturbations are linearly 
additive. Thus these source suppression effects can be addressed in a linear manner, making it possible to 
remove their signature from helioseismic analyses. 

The decreas e of acous tic power in a sunspot has been widely observed and explanations offered (e.g. 
Lites et al.lll982l ): recently. IParchevskv fc Kosovichevi (|2006l ) have suggested that the suppression of convec- 
tion (and hence wave sources) is sufficient to explain more than half of the decrease in acoustic power in 
sunspots. However, in our calculations, even when the source strengths in the region of suppression are 
reduced to zero, we see only about 20% reduction in acoustic power. In any case, it is difficult to compare 
these two results because of the differences in damping rates, the time length of the calculations, the sizes 
of the computational domai ns, etc. To compensate for travel-time measurement 'errors' related to the local 
reduction in acoustic power, iRajaguru et alJ (|2006h proposed a power correction method which we incorpo- 
rated before computing travel times. Since we use broad phase-speed filters and because the acoustic power 



-7- 



is reduced by only 20%, the power correction does little in the way of altering time shifts (~ 10 % change 
at the most) in our simulations. 

Inversions of the mean time shifts (absorbing boundary c ase) using sound-s peed kernels in the ray 
approximation and the multi-channel deconvolution algorithm ( Jensen et al. 19981) are shown in Figure [6] 
The perturbation appears as a shallow (« 1.5 Mm, abutting the photosphere), 7.5 % increase in Sc 2 /c 2 , 
where c is the sound speed. The horizontal size of the anomaly is comparable to that of the region of 
suppressed sources, i.e. about 20 Mm. 
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Fig. 3. — Average cross correlation, C(A, i) for ingoing (on the left) and outgoing (on the right) waves 



obtained from a center-to- annulus scheme ( Duvall et al 



199 



i) for A 



24.35 Mm from simulations with 
absorbing horizontal sides. The solid line shows the average cross correlation for a simulation with a spatially 
homogeneous source distribution ('quiet') and the dashed line for the source-suppressed region. The averaging 
is performed over a region within the 10 Mm disc in the quiet and perturbed simulations. The slight difference 
in amplitudes (there are no phase differences) between the in- and outgoing wave cross correlations in the quiet 
simulation is due to the absence of incoming waves from outside the boundaries. For the source-suppressed 
case, the outgoing wave cross correlation shows a phase advance (roughly 6 seconds) while the corresponding 
ingoing wave correlation shows a much smaller phase s hift. This may contribu te to the asymmetry between 
ingoing and outgoing waves observed in sunspots (e.g., Lindsev fc Braun 20051 ). 



4. DISCUSSION 

We demonstrate that obtaining meaningful travel times is strongly incumbent upon the homogeneity of 
sources in the medium. Numerical and analytical experiments where sources sources were suppressed over a 
region typically the horizontal size of a sunspot predict significant wave phase shifts. Therefore, our analysis 
seems to indicate that helioseismic investigations into the internal constitution of a sunspot are i ncomplete 
witho ut taking into account the effects of inhomogeneously distributed sources and damping jWoodard 



19971 ). We see that in- and outgoing waves are differentially affected, with the asymmetry exacerbated at 
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Fig. 4. — Azimuthally averaged outgoing (solid line) and ingoing (dashed line) time shifts, St, of waves 
that travel distances A = 6.2 Mm (left column) and A = 24.35 Mm (right column). The zero point is the 
center of the source suppression region. The first row (panels a, b) shows time shifts from the simulation 
with absorbing boundaries while the theory of i j2.ll with high damping rates predicts those of the second 
row (panels c, d). The third row (panels e, f) is fro m the simulation with p eriodic boundaries while the 
bottom row (pan els g, h) is from th e theory with the iKorzennik et al.l (|2004r ) damping rates (roughly half 



the linewidths of lDuvall et al.l ll998). A close correspondence is seen between the upper and lower pairs of 



rows. 



increasing travel distance, A, especially when damping rates are high. The large negative mean travel-time 

shifts seen at the shortest travel distances (~ -10 15 s, A = 6.2 Mm) are worrisome for it is not clear 

how accurate estimates are of the amplitude of the near-surface sound-speed perturbation below a sunspot. 
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Fig. 5. — Noise-subtracted (jHanasoge et al.l 120071 ) outgoing (panel a) and ingoing (panel b) time shifts for 
A = 24.35 Mm. We performed three simulations, (I) sound-speed perturbation (of amplitude 7.5% in Sc 2 /c 2 
and size, 20 Mm) + source suppression, (II) only the sound-speed perturbation, and (III) only sources 
suppressed. The perturbations in II and III were identical to the individual components of I. The noise- 
subtracted travel times associated with I (solid line) is seen to be almost indistinguishable from II + III 
(dashed line), indicating that these wave field perturbations are entirely decoupled. 



Similarly, the systematic difference travel times observed for waves (also ~ -6 — -15 s, A = 24.35 Mm) 
that travel larg e r dist ances indicates that the flow inversions may be inaccurate. The power correction of 
Raiaguru et al.l (|2006r ) in this case decreases the magnitude of the travel-time shift s at most b y about 10%. 
The sensitiv ities of other methods of helioseismology like ring diagram analysis ( Hilll X988) and acoustic 
holography (jLindsev fe Braunlll997l ) to the homogeneity of the wave field remain to be investigated. 

The upshot of these calculations is that the methods of Gizon fc Birch (2002) and Birch et al. ( 2004 ) can 
be applied to infer and model out to a large extent the measurement biases introduced by the suppressed 
sources. Moreover, numerical forward models of the solar wave field have beco me increasingly sophisti- 



cated (e.g.. |Parchevskv fe Kosovichevil2006l : ICameron et al.ll2007 : lHanasogelE007bh . presenting ways to test 



inversion results. 



Altered solar model 



Here, we describe the artificially convectively stabilized model (Appendix A of iHanasoge et al.l 120061 ) 
used in our computations. The dimensionless radial co-ordinate is denoted by r, where r expresses fractions 
of the solar radius R@ = 6.959894677 x 10 10 cm. For r < 0.98, background properties as prescribed by model 
S ( Christensen-Dalsgaard et al. 19961) are used. In the range 0.9998 > r > 0.98, the empirical formulae: 

p Q = 4.1522194 [0.998989- r + 4.36138(r-0.98) 21 ] 2 009828 , (Al) 
pa = 2.7392767 x 10 15 [0.998989- r + 4.36138(r- 0.98) 21 ] 3 009828 , (A2) 
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Fig. 6. — Inversion of the noise-subtracted (jHanasoge et al.ll2007fl mean time shifts arising from the sup- 
pressed sources. Shown are piece-wise constant slabs, averaged over the depth range [-0.62, 0] Mm (left 
panel) and [-1.42, -0.62] Mm (right panel), where indicates the surface. Because the inversion is noisy, the 
appearance of random features is observed. Travel-time signatures of suppressed sources and local increases 
in the sound speed look unexpectedly identical, showing significant cross talk from one effect onto the other. 



1 dp 
PoRq dr ' 
Ti = max(lf, 1.507550), 

are implemented, whereas in the region 1.002 > r > 0.9998, an isothermal layer is utilized: 

p = 4.5260638 x 10" 7 exp[7690.7995(0.9998 - r)} 
Po = 1.0252267 x 10 5 exp[7690.7995(0.9998 - r)} 
g = 24998.23 



(A3) 
(A4) 



(A5) 
(A6) 
(A7) 



Density (po) is expressed in units of g cm 3 , pressure (po) in dynes cm 2 , gravity (g) in cm s 2 , the first 
adiabatic index (Ti) is dimensionless, and the sound speed (c) in units of cm s _1 is given by: 



Po 



(A8) 



The eigenfrequencies of the altered model have been computed independently using a boundary value solver 
provided in MATLAB and compared with those recovered from the computations (Figure [T]). The agreement 
is good. 
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B. Code verification 



The accuracy of the numerical scheme described in fJ2] is confirmed using a number of tests (jHanasoge 

2007a)). Before delving into the verification details, it is important to understand the parameter regimes of 
the waves and the limiting factors controlling the simulation timcstcp. The highest frequency of waves of 
interest to us is of the order of 6 mHz, corresponding to a timescale of about 167 seconds. The simulation 
timestep of 2 seconds is significantly smaller than the period of the oscillations. The calculations are evidently 



tempo rally highly over re s olved ; compared to the 4-10 points per wavelength (ppw) quoted by iHu et al 



(|l996l ) and berland et all |2006h . the simulations operate at between 80-250 ppw. In the radial direction, 
the eigenfunctions of the modes contain a rather small number of nodes (10 - 30 depending on the mode) in 
comparison to the actual number of grid points (300 points). The reason for the excessive spatial resolution 
is the need to capture the rapid density (pressure) variation with radius. Therefore, the limiting factor 
controlling the timestep is the large number of density (pressure) scale heights in the computational domain, 
which is why the radial and temporal resolutions are so high. 

We show in Figure [7| that the boundary conditions cause the error convergence rate of the compact finite 
differences to drop to fifth order. Although n ot presented here, the con vergence rate is entirely unchanged 
when the radial de-aliasing filter, described in lHanasoge fc Duvalll (|2007l ). is applied in conjunction with the 
finite differences. Next, to demonstrate the accuracy of the spatial scheme in its entirety (i.e., when used 
with radial de- aliasing and the temporal scheme), we simulate the 1-D propagation of a Gaussian wavelet in 
a box with reflecting boun dary conditions. The g rid-spacing in the calculation follows the constant travel- 
time criterion developed in iHanasoge et ah ( 2006 ). The background model is chosen to be an adiabatically 



stratified, truncated polytrope with index m = 1.5, gravity g = —2.775 x 10 4 cm s 2 e z , reference pressure 



p re f = 1.21 x 10 dynes cm 2 and reference density p re f 
density variations are given by, 



2.78 x 10 7 g cm , such that the pressure and 



p (z) =p re f 



m+l 



and 



Po(z) = p re f ( 

^0 



(Bl) 



(B2) 



The photosphcric level of the background model is at z — 0, with the upper boundary of the truncated 
polytrope pla ced at a depth of zq = 768 km. This model is similar to the stratification in the outer layers of 



the Sun (e.g., iBogdan k Calhlll995h . Because error convergence rates are very sensitive and easily masked 



by slight differences such as the locations of the comparison points of solutions, we start with a highly 
resolved 721 point grid and downsample by successively higher rates (every second point, every third point, 
and so on). The solutions obtained on this sequence of grids are compared with the highly resolved case 
to obtain the error convergence rate. The lower boundary of the simulation is placed at z — —20.876 Mm, 
with wall-like boundary conditions on both ends (v = 0, d z p = —pg, at the boundaries). The timestep of the 
simulation was chosen to be At = 0.05 seconds. The experiment is graphically displayed in Figure[8]and the 
error convergence rate is shown in Figure [9] 
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B.l. Eigenfunctions 

For the polytrop e described above, it is possible to determine the eigenfunctions analytically (e.g., 
Boedan fe Callvlll995h . This will assist us in verifying that the spatial scheme is able to recover the eigen- 



functions accurately. The first step is to set down the equations to be solved: 

d t p(z,t) = -d z (p v) (B3) 
p d t v(z,t) = -d z p-pg (B4) 
d t p(z, t) = -clpod z v + p vg, (B5) 

where p refers to density, c refers to sound speed, the subscript refers to background properties of the 
model, and t time. Differentiating equation (|B4|) with respect to time and substituting for time derivatives 
of density and pressure from equations (|B3[) and (|B5[) respectively, we obtain the following: 

Podfv(z,t) = -d z (-clp d z v + p vg) + d z (p gv). (B6) 

Next we define the Eulerian pressure and velocity fluctuations to be, respectively: 

p(z,t) = p*(z)e-^ (B7) 
v(z,t) = v*{z)e~ luJt . (B8) 



Substituting these expressions into equation (|B6[) . we have: 

- to 2 p z 2 v* = d s ((%p d s v*), (B9) 

where once again, s = —z/zq, pq = p c s m , po = PcS m , Cq = c 2 s, and p c ,p Cl c 2 are the density, pressure 
and sound speed square at s = 1. The upper and lower boundaries of the polytrope are at spatial locations 
s = 1,1?, with D a free parameter. Equation (|B9[) is simplified to obtain: 



2 

sd 2 v* + (m + l)d s v* + ^-v* = 0, (BIO) 

where a — Iujzq /c. Equation (|B10|) is solved to obtain the analytical expression for the eigenfunction: 

v* = As- m/2 J m {as 1/2 ) + Bs-^Y^as 1 / 2 ). (Bll) 

The constants A and B are determined by enforcing the boundary conditions v*(s = 1) = and v*(s — 
D) = 0. From these conditions emerge a sequence of resonant frequencies, a, which can then be used to 
obtain the eigenfunctions of the resonant modes. The eigenfunction for pressure is related to the one for 
velocity according to: 

p* = ^^s m [mv* + sd s v*}. (B12) 
a 

To obtain eigenfunctions from the computations, we excite waves and simulate wave propagation in the 
above-described cavity. Temporal transforms of the entire dataset are computed at each spatial location; 
resonant modes are then isolated by analyzing large amplitude regions in the power spectrum. These 
frequencies are compared to the analytically predicted values to ensure that these are indeed resonant 
modes. Having done so, the temporal spectrum is multiplied by a frequency-window function to retain 
power only in the region of interest and then inverse Fourier transformed. The resultant transforms are the 
desired eigenfunctions. However, spatial error convergence rates are difficult to measure from this experiment 
because the eigenfunction signal is diluted by neighboring modes due to the finite temporal window of the 
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simulations. Moreover the accuracy with which the resonant frequency can be measured is bounded by the 
time length of the calculation. For the eigenfunction shown in Figure [TU1 a resonant mode with v = 6.6111 
mHz was isolated using an extremely narrow, four-point box-car frequency filter. Simulations with varying 
grid spacings all showed a peak in the power spectrum at frequency of 9 /iHz away from the analytical 
prediction (frequency resolution ~ 22/xHz, from a 12-hour simulation). 



B.2. 



As described in Hanasogc ct al 



Efficacy of the transmitting boundary 

()2006h. we use the transmitting boundary conditions of iThompsonl 
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( 1990 ) with an adjoining sponge (e.g., Lui 2003f) to 'prepar e' the waves for the b ound ary. The main r eason 
for using this prescription as opposed to other possibilities (jPoinsot fc Leld (jl992l ); see lColoniusI (]2004r ) for a 
review) is the ease of implementation and efficiency of the method. In the simulations, we use the following: 

(B13) 
(B14) 
(B15) 
(B16) 

with the velocity derivatives computed in a Dirichlet sense, using the values at the end points. 

To test if these boundary conditions change the eigenfunction in any significant manner and to ensure 
that to a large extent, they are indeed non reflecting, we perform ID calculations of wave propagation in a 
background similar to that of § IB. II To give this problem a realistic spin, we stitch an isothermal atmosphere 
to the polytrope so that a finite acoustic cut-off frequency is achieved, thereby providing a natural reflection 
region for the waves. Moreover, we relax the zero-velocity condition on the upper boundary and implement 
a combination of the sponge and transmitting boundary conditions (Eqs. |B13j - |B16j ) while still enforcing 
a zero- velocity condition on the lower boundary. Waves whose frequencies are lower than the acoustic cutoff 
arc reflected back into the interior while an evanescent non-propagating region develops in the isothermal 
atmosphere. Thus, we can determine the effect of the boundary conditions on the simulated eigenfunctions 
by comparing them with their analytical counterparts. 



B.3. Evanescent behavior 

Consider an isothermal layer with constant sound-speed cq with exponentially decaying background 
density (p e )and pressure (p e ) profiles smoothly connected to the truncated polytrope of £|B.1I We have: 

p ref e-^ +z ^ H , (B17) 
Pre f e-( Z0+Z)/H , (B18) 

Tref, (B19) 



Pe = 
Pe = 
T e = 



with z = corresponding to the 'photosphere' of this model, and H to the scale height in the atmosphere. The 
governing equation (|B9[) is unaltered except for the background density and sound speed. Again, we define 
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v(z, t),p(z, t) as: 



p(z,t) 
v(z,t) 



When the constituent equation (|B9[) is solved, we obtain the following for p* and v* 



Pe 

* 



with constants C, D and A a solution of: 

A 



H 



De Xz , 



0, 



1 

2H 

2H 



1-4/1- 



(B20) 
(B21) 

(B22) 
(B23) 



(B24) 

(B25) 
(B26) 



We retrieve two solutions for A and reject the one whose energy density oc pv 2 grows without bound as a 
function z. In this situation, the relation between p* e and v* is given by: 



PcV 
2 i 



V = c A - g 



(B27) 
(B28) 



For boundary conditions, we use normal velocity and Eulerian pressure matching across the boundary 



1: 



v 

_* 
P 



(B29) 
(B30) 



where v* and p* are the velocity and pressure in the polytropic layer, given by equations (|BTT|) and (|B12|) 
respectively. When writing the velocities in the following form, we will have only the pressure equation to 
solve: 



= ^ e -^° s -V 2 [J m («s 1/2 ) + PY m (as^)}, 

PcV 

= A —e~ Xszo {Jm(a)+l3Y m (a)}, 

PcV 



(B31) 
(B32) 



where /3 is the unknown constant we must determine (A can be arbitrarily chosen). Equations (|B27|) 
and (|B32[) constrain p*: 

p* e = Ae- Xsz "+ sz "/ H [J m (a) + /3Y m (a)}. (B33) 



Matching p* = p* at s = 1 gives us the following relations: 

Jm(a) + nJ m -i{a) 



ft 



Y m (a) + KY m -i{a) 
ar\ 



(B34) 
(B35) 
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To determine the resonant modes a of this model, we use the definition of j3 from equation (|B34p and set 
equation (|B3ip to zero at s = D. Having then recovered the resonant frequencies, the pressure and velocity 
eigenfunctions in the interior (s < 1) may be obtained by evaluating: 

v* = A— e- A ^ +2 °/ H S - m / 2 [J m (a S 1 / 2 )+ / 9r m ( as 1 / 2 )], (B36) 

PcC 

p* = -Ane- Xzo+Zo ' H s {m+1]/2 [J m ^ 1 {as 1/ ' 2 ) + l3Y m ^ l {as 1 ' 2 )]. (B37) 



The acoustic-cutoff frequency, lo c , of the model (D > s > 1) is given by: 

CQy/m 2 + 1 1 



(B38) 



The model for this particular test is parametrized by m = 1.5, zq = 768 km, D — 90.6198, cq = 8.51 km 
s _1 , po = 1.21 x 10 5 dynes cm -2 , p = 2.78 x 10" 7 g cm~ 3 , H = z /(m + 1) km, and g = 14160. x 10 5 cm 
s~ 2 . Plotted in Figure fTTI are the analytical (dotted line) and the simulated (solid line) eigenfunctions. The 
sponge is placed adjacent to the upper boundary (located 1232 km above zq). As can be seen the presence 
of the sponge does not affect the interior parts of the acoustic eigenfunction. There is an amplitude error 
near the upper-most region of the polytrope due to the combined influence of the boundary condition and 
the sponge but the nodes remain mostly unaffected. 

A rough test of the efficacy of the boundary conditions is shown in Figure [T2l where an initial Gaussian- 
shaped velocity impulse is allowed to propagate outward. Panel a shows the situation at t = 10 min, and 
the successive panels show the impulses at later instants in time. The amplitude in panel d is of the order 
of 10 -6 , significantly smaller than in panels a through c. Together with the test of Figure QTJ the boundary 
condition seems to allow outward propagating waves to exit the computational domain while leaving the 
eigenfunctions relatively undisturbed. A check of this sort was applied to choose the sponge for the real 
simulations. Since the polytrope + isothermal stratification near the surface is very similar to the model 
used in the simulations, and since the sponges are quite similar in structure, we expect that the eigenfunctions 
in the simulations are also well retrieved while the sponge damps the outward propagating waves. 
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Convergence rate of finite differences 




1.9 2.0 2.1 2.2 2.3 
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Fig. 7. — Spatial convergence rate of the compact finite differences with fifth-order accurate boundary 
conditions. The solid line shows the accuracy of the scheme, while the dashed line is the theoretical fifth- 
order accuracy curve. 
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t = 0, n=721 t = 2 min, n=72l 




0.975 0.980 0.985 0.990 0.995 0.975 0.980 0.9B5 0.990 0.995 



t '. 2 min, ?i 12'., filtered t = 12 min, n 121, unfiltered 




Fig. 8. — Experiment to determine the spatial error convergence rate. The initial condition, a Gaussian 
wavelet in velocity, is shown in panel (a). In (b), the temporally evolved wavelet at time t = 2 min is 
displayed. Simulations are performed with varying numbers of grid points, n — 721, 361, 181, 145, and 121, 
so that each grid is a downsampled version (i.e., every other point, every third point etc.) of the n = 721 
case. Errors are computed at t — 2 min using a downsampled version of the n = 721, t = 2 min solution as a 
template (panel b). In panels (c) and (d), the differences between the n = 121 solution and the downsampled 
n = 721 template at t — 12 min are displayed. The wavelet has propagated all the way out to the upper 
boundary at this point; it is seen that the difference, interpreted as the error, is greater in the unfiltered 
case in panel (d) th an in the filtered version in panel (c), where the filter is applied to dealias variables in 
the radial direction (jHanasoge fc Duvallll2007l ). The difference between (c) and (d), which although appears 
harmless, continues to grow, eventually overwhelming the simulation unless a de-aliasing filter is applied 
frequently. 
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Spatial error convergence rate 
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Fig. 9. — Spatial error convergence rate (with radial dealiasing) based on the experiment of Figure the 
time step was At = 0.05 seconds. The solid line is the error of the compact finite differences and the dashed 
line is a theoretical sixth-order accuracy curve. It is somewhat surprising that the scheme obeys a sixth-order 
accuracy law despite the use of fifth-order boundary conditions. Partly, the reason could be that the problem 
is a consistent initial-boundary value problem, i.e. v = and d z p = —pg at the boundaries. 
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Kio-en function comparison 
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Fig. 10. — Comparison of eigenfunctions for a resonant mode of frequency v = 6.6111 mHz, obtained 
analytically (solid line) and through simulation (dot-dash line) with n = 121. At higher resolutions, the 
two curves are virtually indistinguishable and hence are not shown here. Eigenfunctions for other resonant 
frequencies have also been compared and found to be in good agreement. Including the two boundaries, the 
eigenfunction contains only eleven nodes, far smaller than the number of grid points. With fewer (< 80) 
points, the system develops instabilities because of the steep density gradient. 
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Eigenf unction comparison 
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Fig. 11. — Simulated (solid line) and analytical (dot-dash line) eigenfunctions for v = 1.68 mHz, for the 
model described above. It is seen that the boundary conditions and sponge do not affect the eigenfunction 
over the region of interest; although there is an amplitude error of a few % in the upper-most layers of the 
polytrope, the interior nodes are oblivious to the boundary conditions. This eigenfunction was obtained 
from a 24-hour simulation wherein the waves were constantly excited over a small region in the interior. 
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Fig. 12. — Efficacy of the transmitting boundary. The initial condition is a Gaussian-shaped velocity impulse. 
Panel a shows the situation at t = 10 min, and the successive panels show the impulses at later instants 
in time. The amplitude in panel d is of the order of 1CP 6 , significantly smaller than in panels a through c. 
Together with the test of Figure [Til the boundary seems to do a relatively good job of removing outward 
propagating waves while the interior portion of the eigenfunction is seen to be mostly undisturbed. 



